/*
*		
*		A compilacao e feita por: gcc -o borwein_gmp borwein_gmp.c -lgmp -lm 
*
*/


#include <stdio.h>
#include <time.h>
#include <limits.h>
#include <gmp.h> // gcc -lgmp

#define N 25
#define BASE 10
#define MANTISSA 1024*1024*32


// calcula a_zero e recebe a variavel 'rop' por parametro
// a_zero = 6.0-(4.0*sqrt(2));
void a_zero(mpf_t rop)
{ 
	mpf_t aux;
	mpf_init(aux);
	
	// aux = sqrt(2)
	mpf_sqrt_ui(aux, 2); 
	// aux = aux * 4
	mpf_mul_ui(aux, aux, 4);
	// rop = 6 - aux
	mpf_ui_sub(rop, 6, aux);
	
	// libera a memoria
	mpf_clear(aux);
}


// calcula y_zero e recebe a variavel 'rop' por parametro
// y_zero = sqrt(2)-1.0;
void y_zero(mpf_t rop)
{ 
	mpf_t aux;
	mpf_init(aux);
	
	//aux = sqrt(2)
	mpf_sqrt_ui(aux, 2); 
	// rop = aux - 1
	mpf_sub_ui(rop, aux, 1);

	// libera a memoria
	mpf_clear(aux);
}


// calcula y_kn e recebe as variaveis 'rop' e y por parametro
//	yn = (1-(1-y^4)^1/4)/
//	     (1+(1-y^4)^1/4)
void y_kn(mpf_t rop, mpf_t y)
{ 
	mpf_t y4, n, d;	
	mpf_init(y4);	
	mpf_init(n);
	mpf_init(d);
	
	
	// y4 = y^4
	mpf_pow_ui(y4, y, 4);
	// y4 = 1 - y4
	mpf_ui_sub(y4, 1, y4);
	// y4 = sqrt(y4)
	mpf_sqrt(y4, y4);
	// y4 = sqrt(y4)
	mpf_sqrt(y4, y4);
	
	// n = 1 - y4
	mpf_ui_sub(n, 1, y4);

	// d = 1 + y4
	mpf_add_ui(d, y4, 1);

	// rop = n / d
	mpf_div(rop, n, d);

	// libera a memoria
	mpf_clear(y4);
	mpf_clear(n);
	mpf_clear(d);
}


// calcula a_kn e recebe as variaveis 'rop',a,yn e k por parametro
//an = a(1+ yn)^4 - 2^(2k+3) *  yn(1 + yn + yn^2)
//**	an converge para 1/pi
void a_kn(mpf_t rop, mpf_t a, mpf_t yn, unsigned long k)
{
	mpf_t a4, exp, yk;
	
	
	mpf_init(a4);
	// a4 = yn  1
	mpf_add_ui(a4, yn, 1);
	// a4 = a4^4
	mpf_pow_ui(a4, a4, 4);
	// a4 = a * a4
	mpf_mul(a4, a, a4);

	
	mpf_init(exp);
	// exp = 2
	mpf_set_ui(exp, 2);
	k = (2k+3)
	k *= 2;
	k += 3; 
	// exp = exp^k
	mpf_pow_ui(exp, exp, k);	

	
	mpf_init(yk);
	// yk = yn^2
	mpf_pow_ui(yk, yn, 2);
	// yk = yk + yn
	mpf_add(yk, yk, yn);
	// yk = yk + 1
	mpf_add_ui(yk, yk, 1);
	// yk = yk * yn
	mpf_mul(yk, yk, yn);	
	

	// an = exp * yk
	mpf_mul(exp, exp, yk);
	// rop = a4 - exp
	mpf_sub(rop, a4, exp);
	
	// libera a memoria
	mpf_clear(a4);
	mpf_clear(yk);
	mpf_clear(exp);
}


int main(int argc, char* argv[])
{
	unsigned long int k;
	time_t t_zero, t_final;
	double delta_t;

	// marca o tempo
	time(&t_zero);

	// Float com precisao arbitraria
	mpf_set_default_prec(MANTISSA);
	mpf_t a, y, pi;

	// inicializa as variaveis
	mpf_init(a);
	a_zero(a);

	mpf_init(y);
	y_zero(y);
	
	
	/*

	a = 6 - (4 * 2^1/2)
	y = (2^1/2) - 1

	yn = (1-(1-y^4)^1/4)/
	     (1+(1-y^4)^1/4)

	an = a(1+ yn)^4 - 2^(2k+3) * yn(1 + yn + yn^2)
	**	an converge para 1/pi

	*/
	
	for(k = 0; k < N; k++)
	{
		y_kn(y, y);
		a_kn(a, a, y, k);
    }

	// libera a memoria
	mpf_clear(y);
	
	// an converge para 1/pi
	mpf_ui_div(a, 1, a); 
	
	//mostra o valor de pi
	mpf_out_str(NULL, BASE, 0, a);
	
	// libera a memoria
	mpf_clear(a);
	
	// mostra o tempo gasto
	time(&t_final);
	delta_t = difftime(t_final, t_zero);
	printf("\nTempo = %.2lf segundos\n", delta_t);	

	return 0;
}
